R script to process the raw WTC-3 data to estimate the carbon pools and fluxes for Data Assimilation (DA)

Briefly describe the Carbon Balance Model (CBM)

We planned to use a DA-modelling framework, similar to that used by Richardson et al. (2013) and Mahmud et al. (in preparation). This approach uses a simple carbon balance model shown in Figure 1. The model is driven by daily inputs of gross primary production (GPP). Total maintenance respiration (Rm,tot) is immediately subtracted, and the remainder enters a non-structural C pool (Cn). This pool is utilized for growth at a rate k. Of the utilization flux, a fraction Y is used in growth respiration (Rg), and the remaining fraction (1-Y) is allocated to structural C pools (Cs): among foliage, wood and root (Cs,f, Cs,w, Cs,r). The foliage and root pools are assumed to turn over with rate sf and sr respectively. We assume there is no wood turnover as evidenced from the experiment.

Figure 1: Structure of the Carbon Balance Model (CBM). Pools, shown as solid boxes: Cn, non-structural storage C; Cs,f, structural C in foliage; Cs,r, structural C in roots; Cs,w, structural C in wood. Dummy pools, shown as dashed boxes: Cf,lit, total C in foliage litterfall; Cr,lit, total C in root turnover. Fluxes, denoted by arrows, include: GPP, gross primary production; Rm,tot, maintenance respiration; Rg, growth respiration. Fluxes are governed by seven key parameters: k, storage utilization coefficient; Y, growth respiration fraction; af, allocation to foliage; aw, allocation to wood; \(a_r = (1-a_f-a_w)\), allocation to roots; sf, foliage turnover rate; sr, root turnover rate.

The dynamics of the four carbon pools (Cn, Cs,f, Cs,w, and Cs,r) are described by four difference equations:
\[\Delta C_n = GPP - R_{m,tot} - kC_n\] \[\Delta C_{s,f} = kC_n(1-Y)a_f - s_fC_{s,f}\] \[\Delta C_{s,w} = kC_n(1-Y)a_w\] \[\Delta C_{s,r} = kC_n(1-Y)(1-a_f-a_w) - s_rC_{s,r}\]

Where k is the storage utilization coefficient; Y is the growth respiration fraction; af, aw, ar are the allocation to foliage, wood and root respectively; sf and sr are the leaf and root turnover rates respectively; \(\sum s_fC_{s,f} = C_{f,lit}\) and \(\sum s_rC_{s,r} = C_{r,lit}\) are the dummy foliage and root litter pools respectively.

Total maintenance respiration, Rm,tot was calculated as a temperature-dependent respiration rates for foliage, wood and root (Rm,f, Rm,w and Rm,r respectively), multiplied by plant organ C masses (Ct,f, Ct,w, and Ct,r are the total C in foliage, wood and root respectively). Wood respiration was further partitioned into stem and branches. Similarly root respiration was partitioned into different root size classes (fine, intermediate, coarse and bole roots). Growth respiration (Rg) at each time step was calculated as a modeled respiration rate (Y) multiplied by plant biomass changes at that time step (\(\Delta C_{t,f} + \Delta C_{t,w} + \Delta C_{t,r}\)).

\[R_{m,tot} = R_{m,f} C_{t,f} + R_{m,w} C_{t,w} + R_{m,r}C_{t,r}\] \[R_{m,w} C_{t,w} = R_{m,s} C_{t,s} + R_{m,b} C_{t,b}\] \[R_{m,r} C_{t,r} = R_{m,fr} C_{t,fr} + R_{m,ir} C_{t,ir} + R_{m,cr} C_{t,cr} + R_{m,br} C_{t,br}\] \[R_g = Y(\Delta C_{t,f} + \Delta C_{t,w} + \Delta C_{t,r})\]

Where Rm,s, Rm,b, Rm,fr, Rm,ir, Rm,cr and Rm,br are the maintenance respiration rates with Ct,s, Ct,b, Ct,fr, Ct,ir, Ct,cr and Ct,br are the total C in stem, branches, fine roots, intermediate roots, coarse roots and bole roots respectively.

The non-structural (storage) C pool (Cn) is assumed to be divided among foliage, wood and root tissues (Cn,f, Cn,w, Cn,r) according to empirically-determined fractions.
- However, WTC-3 experiment only measured leaf non-structural C (Cn,f), and therefore to estimate the partitioning of the non-structural C among different organs, we hope to use data from WTC-4 experiment on similar-sized seedlings of a related species (Eucalyptus parramattensis). We will consider different treatments from the experimental dataset, to find the Cn partitioning to foliage, wood and roots.
- Another possibility (if WTC-4 data are not available!!) would be to assume another two parameters to estimate the Cn partitioning to foliage, wood and roots (pf, Cn partitioning to foliage; pw, Cn partitioning to wood; \(p_r = (1-p_f-p_w)\), Cn partitioning to roots).

Total carbon in each tissue (Ct) is then calculated as the sum of non-structural carbon (Cn) and structural carbon (Cs) for that tissue.

\[C_{t,f} = C_{n,f} + C_{s,f}\] \[C_{t,w} = C_{n,w} + C_{s,w}\] \[C_{t,r} = C_{n,r} + C_{s,r}\] We plan to estimate seven parameters (k, Y, af, aw, ar, sf, sr) of the CBM for this experiment using DA. GPP and maintenance respiration rates will be used as model inputs in the DA framework, whereas the measurements of total aboveground respiration (Rabove), total C masses (Ct,f, Ct,w, Ct,r), foliage litterfall (Cf,lit) and foliage NSC (Cn,f) will serve to constrain the parameters. Aboveground respiration (Rabove) was estimated by summing both maintenance and growth respiration components of foliage and wood.
\[R_{above} = R_{m,f} C_{t,f} + R_{m,w} C_{t,w} + Y\Delta C_{t,f} + Y\Delta C_{t,w}\]

Step 1:

Select treatment strategy
  • Check whether there are any pre-existing differences between the chambers of droughted treatments
  • Should we consider:
    1. the droughted treatments separately from the start of the flux measurements (Sept 2013): n = 3 or
    2. from actual drought implementation on 12 Feb 2014: n = 6 predrought and n = 3 postdrought ??
Data used:
  1. WTC_TEMP_CM_TREE-HEIGHT-DIAMETER_20121120-20140527_L1_V2.CSV

Figure 2: Growth (Diameter) of Eucalyptus tereticornis trees exposed to warming and drought

Finding:

It is obvious from Figure 1 that there were pre-existing differences between the warmed watered (black vs blue, Figure 1) and warmed drought (black vs pink, Figure 1) trees. It seems like the warmed-drought trees actually started the drought smaller than the warmed-watered trees (blue vs pink, Figure 1). If we go for the above option 2 (i.e. separate all 4 treatments from actual drought implementation on 12 Feb 2014), the diameter time series (red line, Figure 1) makes it look like there was a strong effect of the drought on the tree growth but it was more a function of pre-existing differences in tree size. The height data showed the similar pattern.

So based on both diameter and height data, we decided to represent drought and watered treatments separately from the beginning of the experiment. So total number of treatments = 4 (ambient+watered, ambient+droughted, elevated+watered, elevated+droughted) with sample size of n = 3 for each treatment.

Step 2:

Inputs:
c1 = 0.48 # (unit conversion from gDM to gC: 1 gDM = 0.48 gC)

Step 3: Read and process GPP and Rabove (aboveground respiration) data

Partitioning CO2 fluxes into GPP and Rabove

Partitioning of the hourly net CO2 fluxes into the components of GPP and Rabove were done using a technique common to eddy-covariance research (Reichstein et al. 2005); described thoroughly in Drake et al. (2016, 2017-in preparation). We assumed GPP to be zero at night when PPFD = 0, indicating the measured net C flux in such conditions was used as the measure of Rabove. We utilized the direct measurements of whole-canopy Rabove and its temperature dependence at night to predict Rabove for each hourly measurement as a function of air temperature. We then calculated GPP as the sum of the measured net CO2 flux and the predicted Rabove, given the measured air temperature.

Data used:
  1. WTC_TEMP_CM_WTCFLUX_20130914-20140526_L2_V2.csv

Figure 3: Daily GPP and aboveground respiration (Rabove) of Eucalyptus tereticornis trees for all 4 treatments. The grey bars represent standard errors for sample size, n=3.

Step 4: Estimate biomass pools for each treatment

  1. Process stem height and diameter for various treatment cases
  • Tree height (cm) and stem diameter (mm) were measured fortnightly. Stem diameter was measured at 30-cm-intervals along each tree stem from a basal height of 15-cm (prior to floor installation) or 65-cm (after floor installation) to the tree apex. However, the diameters at 15-cm height are essential to estimate the initial root mass (at the start of the flux measurement), based on the allometry and harvest root mass of Court’s experiment. Therefore, we gap-filled the 15-cm diameters by fitting a linear regression between diameters at 15 and 65 cm heights.
  • Finally we averaged both the height and diameter (at 15-cm) across the treatments (n=3) to get the mean fortnightly allometry data with standard errors.
  1. Estimate the root mass for each treatment
  • We estimated initial root mass (g C) using allometry (height and diameter at 15-cm height) and harvest root mass of control (free) seedlings from Court’s experiment on similar seedlings (Euc. Teri.). We read and merged all necessary data from Court’s experiment to get the allometry and harvest root mass. We considered both data from Court’s experiment and WTC-3 harvest for the control (ambient and watered) treatments, and fitted a linear regression (with interaction) to estimate root mass (with standard errors) from the initial data of stem diameter and tree height at the start of the flux measurements: lm(log10(rootmass) ~ log10(diameter) * log10(height), P < 0.001, r2 = 0.99.
  • We averaged WTC-3 final harvest root mass across the treatments (n=3) to get the mean and standard errors.
  1. Get an estimate of stem, branch and leaf mass for each treatment following the script of John for WTC-3 flux partitioning paper (in progress)
  • Fortnightly stem mass (g C) was calculated from geometry using stem diameter at different heights. We considered stem as i) frustum of a cone above 65 cm, and ii) cylinder below 65 cm to ground.
  • Branch mass was estimated using harvest branch mass, branch census and stem volume. Using the harvest data, we setup a log-log regression for branch mass and diameter: (log10(branch mass) = -1.299 + 2.722 × log10(branch diameter), P < 0.001, r2 = 0.91, branch mass in g, branch diameter in mm, n = 48 branches).
  • Then we applied this allometry-mass regression to the branch census dataset to estimate the branch mass on three intermediate dates (24 Oct 2013, 15 Jan 2014, and 22 May 2014). Total branch mass and stem volume were strongly correlated in a chamber specific manner (log-log ANCOVA, P < 0.001, r2 = 0.95), which was used to estimate branch mass (g C) as a function of stem volume for each fortnightly growth interval.
  • We estimated fortnightly leaf mass using canopy-weighted SLA (Specific leaf area) from harvest and leaf area (LA) estimates according to Barton et al. (2012) and Drake et al. (2016). LA was based on stem height and leaf litterfall at any time. Standing LA was measured twice on 9 Sept 2013 and 10 Feb 2014 for each tree by counting all the leaves and multiplying by a tree-specific mean leaf size measured across the canopy of each tree with a handheld leaf area meter. A third direct measurement of standing canopy leaf area was calculated from the final harvest data (26 May 2014) by multiplying total canopy leaf dry mass by SLA weighted by the leaf dry mass in each canopy layer.
Assumptions:
  1. No LMA variation over time: We plotted the LMA variation over time to see whether there was any time effect on LMA. We used all available data from Mike, Court and harvest to investigate the trend, however did not find any concrete relationship. Therefore, we ended up with the constant LMA from harvest data as used by John and suggested by Remko, Court.

Figure 4: Estimated biomass pool of WTC-3 Eucalyptus tereticornis trees for all 4 treatments. We predicted total C mass present in both foliage and wood with fortnightly interval, however only the initial and final C mass for roots. Note that, the points are jittered in all figures to see the grey standard error bars (n=3).

Step 5: Process the leaf litter data

Litterfall was collected, dried, and weighed fortnightly for each tree. The dummy foliage litter pool, Cf,lit (g C) was estimated by cumulative summing of fortnightly litter starting from the flux measurements. The daily foliage litterfall was predicted as a linear function of time between two fortnightly consecutive measurements and will be used to constrain the parameter foliage turnover rate, sf.

Data used:
  1. WTC_TEMP_CM_LEAFLITTER_20130913-20140528_L1.csv

Figure 5: Dummy leaf litterfall pool of WTC-3 Eucalyptus tereticornis trees for all 4 treatments. Note that, the points are jittered to see the grey standard error bars (n=3).

Step 6: Estimate the leaf TNC (storage) pool

  1. First we considered the foliage TNC measurements over time and averaged the data across the treatments (n=3) to get the mean storage with standard errors.
  2. We then merged the diurnal foliage TNC data by averaging all 2-hourly interval data on a sunny (20 February) and an overcast/rainy day (26 March). Similarly this data were averaged across the treatments (n=3) to get the mean with standard errors.
  3. We also took onto board the diurnal leaf TNC data when trees were girdled on 15 May 2014 (only the ambient watered treatments).
  4. The unit was converted to gC in TNC per gC of plant for consistency.
Data used:
  1. WTC_TEMP_CM_LEAFCARB_20130515-20140402_L2.csv
  2. WTC_TEMP_CM_LEAFCARB-DIURNAL_20140220-20140326_R.csv
  3. WTC_TEMP_CM_PETIOLEGIRDLE-LEAFMASS-AREA-CARB_20140507-20140512_L2.csv

Figure 6: Foliage TNC pool, Cn,f over time of WTC-3 Eucalyptus tereticornis trees for all 4 treatments. Note that, the points are jittered to see the grey standard error bars (n=3).

Step 7: Estimate the partitioning of woodmass to branch/stem and rootmass to fine/intermediate/coarse/bole roots

Purpose: To predict total wood and root maintenance respiration \(R_{m,w}C_{t,w}\) and \(R_{m,r}C_{t,r}\), we need both daily respiration rate and corresponding daily C mass of all wood and root components.
Woodmass partitioning
  1. Calculate partitioning of WTC-3 woodmass: We did not find any significant difference across the treatments for wood partitioning, so we considered equivalent wood partitioning for all treatments (Figure 7, solid lines).
  2. We then considered a linear variation over time in between the fortnightly measurements to predict the daily woodmass partitioning.
Rootmass partitioning
  1. Calculate partitioning of WTC-3 harvest rootmass: We did not find any significant difference across the treatments for root partitioning, so we considered equivalent root partitioning for all treatments (Figure 7, dotted lines end points).
  2. Estimate the partitioning of initial rootmass based on Court’s Pot experiment harvest rootmass (control free seedling): We assume there were only fine and intermediate roots, with no coarse and bole roots at the start of flux measurements (Figure 7, starting points).
  3. We then considered a linear variation over time to predict the daily rootmass partitioning, in agreement with above ground biomass partitioning that follows linear trend.
Data used:
  1. WTC_TEMP_CM_LEAFCARB_20130515-20140402_L2.csv
  2. Below-ground Sink limited container experiment (Court’s exp) - Harvest rootmass data of free seedlings

Figure 7: Rootmass partitioning of WTC-3 Eucalyptus tereticornis trees for all 4 treatments. The grey shade shows the standard error (n=12).

Step 8: Plot leaf-scale respiration vs. leaf temperature to find Q10 for WTC-3

  1. We plotted the leaf-scale respiration (Rleaf) against leaf temperature (Tleaf) to test the short-term temperature sensitivity of the respiration of individual leaves similar to Drake et al. (2016). However, instead of using Arrhenius function (Drake et al. 2016), we fitted the exponential relationship: \(R_{leaf} = R_{leaf,25}Q_{10}^{(T_{leaf}-25)/10}\). The statistical test did not show any significant differences for drought treatment, however warming reduced the pre-exponential coefficient (Rleaf,25) but did not alter Q10, reflecting a Q10 value of 2.26 at 25°C for all treatments that agrees with previous study of Drake et al. (2016).
Data used:
  1. WTC_TEMP_CM_GX-RdarkVsT_20140207-20140423_L1.csv

Figure 8: The short-term temperature sensitivity of the respiration of individual leaves relative to leaf temperature. Error bars reflect the standard errors of Eucalyptus tereticornis trees exposed to each temperature treatments (n = 6).

Step 9: Relationship between mass-based leaf-scale respiration at 25°C and air temperature

  1. We plotted the mass-based leaf-scale respiration at 25°C, Rleaf,25 (nmol CO2 g-1 s-1) over time following Aspinwall et al. (2016). Figure 9 (first panel) shows that leaf R (mass-basis) measured at 25°C varied over time and decreased as mean air temperature (Tair) increased (rest of the panels).
  2. Using the site weather data, we fitted a linear regression between Rleaf25 and Tair. We tried different options with mean air temperature of previous 7-days, previous 3-days, previous day, 3-day including the day of respiration measurement and only the day of respiration measurement. Treatment had no effect on the intercept or slope of the relationship between Rleaf25 and Tair. Based on the data fit and considering the finding of previous study (Aspinwall et al. 2016), we found the best relationship of 25°C mass-based leaf-scale respiration (nmol CO2 g-1 s-1) with 3-day mean air temperature, Tair,3day (°C) with P < 0.001, r2 = 0.75:
    \[R_{leaf,25} = 21.5 - 0.565 T_{air,3day}\]
Data used:
  1. WTC_TEMP_CM_GX-Rdark25_20130617-20140402_L2.csv
  2. WTC_TEMP_CM_WTCMET_20130601-20130630_L1_v1.csv - WTC_TEMP_CM_WTCMET_20140501-20140531_L1_v1.csv

Figure 9: Mass-based night-time leaf-scale dark respiration measured in Eucalyptus tereticornis at a set temperature of 25°C through time (first panel) and in relation to prevailing mean air temperature, Tair (rest of the panels). Mean values are those of replicate whole-tree chambers (n = 6 for ambient/warmed treatment), determined based on three individual leaf measures per chamber. No statistically significant variation for drought/watered treatment on leaf respiration.

Step 10: Estimate daily mean respiration rates for foliage, wood and roots

  1. Foliage respiration: First we calculated the daily mean mass-based leaf-scale respiration, Rleaf,25 (nmol CO2 g-1 s-1) at 25°C using the regression we fitted with 3-day mean air temperature, Tair,3day (°C) in Step 9: \(R_{leaf,25} = 21.5 - 0.565 T_{air,3day}\).
  2. Bole wood respiration: There was no statistically significant difference in bole/stem wood respiration across the treatments, and the rate was 0.656 nmol CO2 g-1 s-1 at 15°C.
  3. Branch wood respiration: We found statistically significant difference in branch wood respiration across the temperature treatments (ambient vs. warmed). The rates of branch wood respiration in WTC-3 were about 1.33 nmol CO2 g-1 s-1 in the ambient treatment and about 1.04 nmol CO2 g-1 s-1 in the warmed (+3°C) treatment when measured at 15°C. This indicates temperature acclimation in branch wood respiration.
  4. Fine root respiration: We assumed the specific rate of fine root respiration at 25°C was about 10 nmol CO2 g-1 s-1 (Drake et al. 2017). We estimated fine root respiration at 15°C using the Q10 value of 2.26 (as found in Step 8).
  5. Coarse root respiration: For coarse roots (Diameter > 10 mm), we suggested using the rate we determined for branch wood in E. tereticornis. The rationale would be that the woody coarse roots were made up of more or less well developed xylem (and phloem) tissues, were thus more or less functioning as underground woody branches. Therefore, we would expect similar temperature acclimation in coarse root as in branch wood respiration.
  6. Bole root respiration: The bole and big tap roots were assigned the bole wood respiration, collected at the end of the WTC-3 study at 15°C.

  7. Intermediate root respiration: For the intermediate size class of 2 to 10 mm, We interpolated respiration rates between the size classes to estimate the middle class using a simple log-log plot.

  8. 15-mins interval respiration rates:
  • We estimated 15-mins interval foliage respiration rates, Rleaf at 15-mins air temperatures (Tair) using the site weather data, mass-based daily mean leaf-scale respiration at 25°C, Rleaf,25 and Q10 of 2.26: \(R_{leaf} = R_{leaf,25}Q_{10}^{(T_{leaf}-25)/10}\). We assume \(T_{leaf} = T_{air}\).
  • Wood respiration rates in 15-mins interval were calculated at air temperature (Tair) using the site weather data, mass-based wood respiration rates at 15°C and Q10 of 2.26.
  • Similarly root respiration rates in 15-mins interval were calculated using soil temperature at 10 cm depth, mass-based root respiration rates at 15°C and Q10 of 2.26.
  1. Daily mean respiration rates: Daily mean respiration rates for all tree components were calculated by summing all 15-mins data for each day.
  2. The units were converted to gC g-1C d-1 (from nmol CO2 g-1 s-1) for consistency.
Data used:
  1. WTC_TEMP_CM_GX-RBRANCH_20140513-20140522_L1_v1.csv
  2. WTC_TEMP_CM_WTCFLUX-STEM_20140528_L1_v1.csv
  3. WTC_TEMP_CM_WTCMET_20130601-20130630_L1_v1.csv - WTC_TEMP_CM_WTCMET_20140501-20140531_L1_v1.csv

Figure 10: Daily mean respiration rates for all tree components (foliage, wood and roots). There were no statistically significant difference in respiration rates across drought/watered treatments. Both branch wood and coarse root respiration showed temperature acclimation having higher rates in ambient than warmed condition.

Step 11: Apply Data Assimilation (DA) with the estimates of WTC-3 carbon pools and fluxes to predict all seven parameters (k, Y, af, aw, ar, sf, sr), with the aim to quantify the main C balance processes (respiration, carbohydrate utilisation, allocation and turnover) in response to elevated temperature and drought on Eucalyptus tereticornis trees grown in WTCs.

Hypothesis 1: NSC vs. Temperature - As Paul, Driscoll & Lawlor (1991) have shown, low T generally results in accumulation of carbohydrates, indicating that growth or storage are limiting. Conversely, at high temperatures, sink demand is large and assimilate is depleted so there should be a marked higher utilization rate (k).
- The carbohydrate concentration decreased in stem and root tissues for Citrus plants, while it increased in leaf tissues under moderate warm conditions (30/20ºC than at 25/20ºC, Ribeiro et al. 2012).
- We hypothesize that drought is carbon limiting and negatively impacts plant carbon balance and that plants will rely on stored carbon to survive carbon limitation and therefore deplete their stored carbon reserves over time.

… in progress …